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Abstract 



We consider the alternating method of Kozlov and Maz'ya for solving the Cauchy 
n ' problem for elliptic boundary-value problems. Considering the case of the Laplacian, 

we show that this method can be recast as a form of Landweber iteration. In addition 
to conceptual advantages, this observation leads to some practical improvements. We 
show how to accelerate Kozlov-Maz'ya iteration using the conjugate gradient algo- 
rithm, and we show how to modify the method to obtain a more practical stopping 
criterion. 



> 

^ 1 Introduction 

(N 

[^ The Cauchy problem for elliptic equations, where Dirichlet and Neumann data are pre- 

O scribed simultaneously on a strict subset of the domain boundary, is a prototypical ill-posed 

^ problem. It is a linear problem, and can be approached using any of a number of standard 

r| regularization techniques such as Tikhonov regularization [CDJPOl], as well as logarithmic 

convexity methods[Pa87]. Kozlov and Maz'ya [KM90] (see also [KMF91]) introduced a 
novel method for solving this problem that, while related to the general class of iterative 
C^ methods, was not evidently one of the standard ones. The method (sometimes known as 

the alternating method and called here Kozlov-Maz'ya iteration) is straightforward to im- 
plement numerically, and is therefore an attractive choice for practical use. There are some 
drawbacks, however. The formal stopping criterion for Kozlov-Maz'ya iteration involves 
error estimates in certain Sobolov spaces with fractional derivatives. Arriving at such es- 
timates for real data poses some difficulty. Moreover, it has been observed that Kozlov- 
Maz'ya iteration suffers from being slow. Although there have been efforts to accelerate 
the method using certain relaxation factors [JN99] [JLM04], formal proofs of the stability 
of these ad-hoc techniques are not available. 

Our primary result is a demonstration that Kozlov-Maz'ya iteration is, in fact, a form of 
Landweber iteration [La51] between function spaces equipped with suitable norms. This 
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observation yields a number of advantages. First, the extensive body of literature con- 
cerning Landweber iteration can be brought to apply to Kozlov-Maz'ya iteration. Proofs 
concerning its stability and rate of convergence can then be quoted from textbooks. Second, 
standard techniques for accelerating Landweber iteration can be applied to Kozlov-Maz'ya 
iteration. We indicate here a variation based on the conjugate gradient method that is nearly 
as simple as standard Kozlov-Maz'ya iteration and leads to a fast, order-optimal regular- 
ization method. Finally, we show how to modify one of the function spaces involved in 
the Landweber and conjugate gradient iterations to obtain a similar method with a more 
practical stopping criterion. 

The motivation for this work comes from an inverse problem in glaciology [MTAS08], 
which considered the Cauchy problem for a nonlinear elliptic PDF. In that paper, the lin- 
earized inverse problems were solved using Kozlov-Mazya iteration, accelerated by the 
techniques described below in Section 5. For simplicity, we focus our attention here on a 
model elliptic problem for the Laplacian; the extension to more general elliptic operators 
is straightforward. We remark that Kozlov-Maz'ya iteration has been extended to certain 
parabolic and hyperbolic inverse problems [BLOl] as well as to a degenerate elliptic prob- 
lem (the Stokes system) [BJKL05]. We do not treat these problems, but it hoped that the 
ideas presented here might also be useful in these cases. 



1.1 Formulation of the model Cauchy problem 

Let n C M" be an open, bounded, and connected set with a Lipschitz boundary dVl. Sup- 
pose S and B are nonempty open subsets of (90 sharing a common boundary 11 and that 
dVl = SUBUH is a Lipschitz dissection as defined in [McOO] (effectively 11 is an embedded 
Lipschitz hypersurface of F). We suppose that boundary data is known on S but unknown 
on B; the notation suggests that S is an accessible surface and B is an inaccessible base. 
The Cauchy problem for the Laplacian is the following: 

— Au = f in r2 

u = a on S (1) 

dnU = T on S. 

Here / G {H-\n))*, a G H^/\S), r G H-^/\S), and dn denotes the the normal deriva- 
tive. 

Our notation and conventions for Sobolov spaces follow [McOO], except for one case noted 
below. The space H'^(S) is the set of restrictions of distributions in H^{dil) to S and has 
the quotient norm; in this paper we will only use the cases s = —1/2, 0, 1/2. The subset 
of distributions a G H'^(d^) such that a\B = is denoted by Hqq(S) (and by H^{S) in 
[McOO]). It is a closed subspace of H'^^dVl) and inherits the norm from the larger space. 
We will consider elements of Hqq(S) as elements of H^{S) or as elements of H'^{dVl) 
interchangeably and without comment. Because of the regularity of the sets S and B, 
there is a natural identification of H^^/'^{S) with the dual space of H^q (S). More details 
concerning these conventions can be found in the Appendix. 



The key step of recasting Kozlov-Maz'ya iteration as Landweber iteration requires making 
a judicious choice of (equivalent) norms on these boundary Sobolov spaces so that the ad- 
joints of certain operators take on a natural form. Proofs of the equivalence of these norms 
can be found in the Appendix. For brevity we describe these norms here for distributions 
on S with obvious adjustments needed for distributions on B. 

1.1.1 The space H-^/\S) 

Let ip e H^^/^{S) and let v^ be the solution of 

A t>^ = in r2 
dnV^ = ip on S 
v^ = on B. 

Then 

Similarly, ifipi,'ip2 G H-^/'^{S) then 



Jn 
This norm was described in the original work on Kozlov-Maz'ya iteration [KM90]. 

1.1.2 The space Hol\S) 

Let (j) e Hll^'^iS), so G H^/'^{dQ) and (f)\B = 0. Let w^ be the solution of 

Aw^ = in fi 
w^ = (j) on dVt. 

Then 

Ii0ii^i/2(s) = / |vw</,|^ 

and there is a corresponding inner product 

1.1.3 The space H^''^{S) 

Let e H^/'^{S). Let r^ be the solution of 

A r^ = in fi 
r^ = (f) on S 
dnf^ = on 5. 



Then 
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where 15*1 = J^l. There is a corresponding inner product defined analogously to the other 
spaces. 



2 Kozlov-Maz'ya iteration 

Kozlov-Maz'ya iteration proceeds by alternating between solving boundary-value prob- 
lems involving the Dirichlet data (a) and Neumann data (r) in turn. Given ip E H^^^'^{B), 
let J\f{ij) = V where v E H^{n) is the solution of 

-Av = f in 17 
V = a on S 
dnV = Ip on B. 

Given E H^/'^{B), let P(0) = w where w E H^{Vt) is the solution of 

— Aw = f in n 
dnW = T on S 
w = (p on B. 

The solvability of these equations is discussed briefly in the Appendix. 

Starting with an initial element ^o £ H^^^'^{B), we let vq = A/'(^o) and wq = V{vo\b). 
Then tpi = 9„Wo|b- Sequences {^/t}, {vk} and {wk} are obtained by repeating these 
operations. Letting 

KM{tlj)=dnV{Af{^)\B)\B, 

we see that ipk+i = KM{il>k). Note that we use the convention that operators with a script 
font yield distributions on Vt, whereas operators with a roman font yield distributions on a 
subset of dVt. 

It was proved in [KMF91] that u E H^{Vl) is a solution of the Cauchy problem (1) if and 
only if dnu\B is a fixed point of KM. Moreover, if there exists a solution u of equation 
(1), then the functions f„ and w„ converge to u in H^iVt). And finally, for approximate 
boundary data {as.rs), stopping the iteration early according to a discrepancy principle 
leads to a regularization strategy for solving the Cauchy problem. 

There is a dual formulation of Kozlov-Maz'ya iteration obtained via the operator MK : 

H^'^{S) ^ H^I^{S) defined by 



MK{<j)) = U{dnV 



B B 



(i.e. MK is KM performed in the reverse order). In the following two sections we will 
show that the original form of Kozlov-Maz'ya iteration and its less- studied dual formula- 
tion can be exhibited as forms of Landweber iteration. 



The maps V, J\f, KM, and MK are affine, and it will be useful to have notation for their 
linear parts. Let Vq, Aq, KMq, and MKq be defined as before but with homogeneous data 

(/ = 0, a = 0, r = 0). It then follows that for any ^, V^ G H-^''^{B) and 0, G H^I'^{B) 

KM{'4) + i)) = KM{ip) + KMo{i)) 
MK{4> + 0) = MK{^) + MKo{(p). 



3 Landweber iteration (Neumann version) 

Define N : H-^/^{B) -^ H-'^/^{S) by 

N{i;)=dnUms. 

We will show in this section that the previously described alternating technique is, in fact, 
the Landweber method applied to the operator equation 

N{^) = T. (3) 

Note that if m is a solution of the Cauchy problem (1), then N{dnu\B) = t. Moreover, if z^ 
solves equation (3), then u = N'i'ip) solves the Cauchy problem (1). 

Let Nq be defined similarly to N using A/q in place of A/" (i.e using homogeneous data). So 
No : H-^/\B) -^ H-^/^{S) is a linear map and 

iV(V') = iV(0) + iVo(^). 

Hence the operator equation (3) can be rewritten 

iVo(^) = r-iV(0). (4) 

The Landweber method provides a regularization technique for solving equation (4) that 
proceeds by minimizing the functional 

using a steepest descent algorithm. The gradient of J dXip is 

-n; [t - iv(o) - N^m , 

and we define 

LnW = V^ + aN* [t - iV(0) - NoW] (5) 

where a is a fixed constant chosen so that < a < 1/| | A^ol P- The Landweber method then 
produces iterates ^Jk+i = Lj^fi^xpk) starting from an initial value ^q, and the functions 4'k 
are then approximate solutions of the original operator equation (3). 

Computation using the Landweber method requires knowledge of the adjoint A^q, which 
has a natural form given our chosen inner products. 



Lemma 1. Let ^ E H ^^'^(8), and define q to be the solution of 

— Aq = in Q 
dnq = C, on S 
q = on B. 

Then N*{0 = dnq\B. 

Proof. Let tp E H^^^'^(B) and ^ E H^^^'^{S) be arbitrary. We then let v, w, q and r be the 
solutions of the following boundary value problems: 

-Ai; = -Aw = -Aq = -Ar = inf] 

V = dnW = dnV dnq = ^ r = on S 

dnV = ip w = g = dnT = dnq on B. 

Then 



(A^oV',Oii'-i/2(5) = / Vw-Vq= / dnwq= / dnW q 
Jn Jdo, Js 



dnV Q = dnV Q = Vf ■ Vg. 
'5 Jdn Jn 

(6) 

The first equality follows from the definition of the inner product on H^^^'^(S) and we have 

used q = on B and dnV = dnW on S in the subsequent equalities. Similarly, 



{lp,dnq\B)H~^/HB)= VvVr = V dnT = V dnT 

Jn Jdn Js 



V dnq = V dnq = / Vf ■ Vg. 
Is Jdn Jn 

(7) 

From equalities (6) and (7) we conclude 

(^0^,0/^-1/2(5) = (V',5ng|ij)H-i/2(B) ' 

and since tp E H^^^'^(B) is arbitrary, Nq{^) = q\B as claimed. D 

The alternating method of [KM90] is exactly the Landweber method (with the constant 
a = 1) applied to minimizing the functional J. 

Proposition 1. For any ip E H^^^^{B), we have 

L^i^) = KM{tP). 

Consequently, the iterates produced by the (Neumann) Landweber method and the (Neu- 
mann) Kozlov-Maz'ya alternating method are identical. 



Proof. Lemmas 2 and 3 proved below establish that 

N;N,{i,) =^- KMoi^) (8) 



and 
Hence 



N*{t-N{0)) = KM{0). (9) 



L^i^) = V^ + iV*(r - iV(0)) - N*Noi^P) 
= ip + KM{0) -[ip- KMoii;)] 
= KM{0) + KMo{^) 
= KM{^) 

where in the last step we have used the decomposition of the affine map KM. D 

It remains to establish equations (8) and (9), which is done in following two lemmas. 
Lemma 2. For any ip in H^^/'^{B), 

N;No{^)=ij-KMo{'ilj). 

Proof. Let ip E H^^^'^{B) and let v and w solve the boundary-value problems 

f = dnW = on 5 

dnV = ip w = V on B. 

Then A^o(^) = dnv\s and KMq{iP) = 9„w|b- On the other hand, let g = f — w. Then 

- A g = in f] 
dnQ = dnV on S 
q = on B 

and hence A* (5„f 1 5) = dnqls- But 

dnq\B = dnV\B " i9„w|b = tp - KMq{iI)). 



So 



Lemma 3. 



iV*iVo(V>)=V'-irMo(^). 



Ar*(r-Ar(0)) =irM(0) 



D 



Proof. Let v and w be the solutions of 

-At; = / -/\w = f mVL 
V = a dnW = T on S 

dnV = w = V on B. 

Then A^(0) = dnv\s and KM(0) = dnw\B- Let r = w — v. Then 

- A r = in O 

dnV = T - dnv\s on S* 
r = on 5. 

By Lemma 1 , 

A^o(^ - dnV\s) = dnr\B = dnW\B - dnV\B- 

But dnv\s = N{0), dnw\B = KM{0) and dnvls = 0. Hence N*{t - iV(0)) = KM{0). 

D 

Remark 1. It follows from Lemma 2 that the operator KMq : H~^/'^{B) — > H^^/'^{B) is 
self adjoint. This fact was proved independently in [KM90], and it played a central role in 
their results. The new observation in the current work is that this self-adjointedness arises 
because of equation (8) and that the Kozlov-Maz'ya iteration is intimately connected with 
a minimization procedure. 

We have now proved all of the ingredients of Proposition 1 . 

To justify the use of Landweber iteration with relaxation factor a = 1 in equation (5), we 
require an estimate of the norm of A^q- 

Lemma 4. The operator norm of Nq satisfies \\Nq\\ < 1. 

Proof. Let v,w, and q be solutions of the following boundary value problems: 

-At; = -Aw = -Aq = infi 

V = dnW = dnQ = dnV on S 

dnV = ip w = V q = on B. 



Then iVo^ = dnv\s = dnq\s and 



H-y^B) = I iVt;P 



||iVoV^||^-i/.(5)= / |Vg|^ 
Jn 

Now let r = f — w. Then r is harmonic, r\B = 0, and 9„r = dnV on S. Hence v — w 
r = q. Moreover, 



/ Vf ■ Vw = V (dnw) = w (dnw) = / |Vw 
Jn Jan Jan Jn 
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since v = w on B and d„w = on S. Hence 



\No^\\l-i/2,s) = /|Vg|^= f\Viv-w)f= [ \Vvf-2Vv-Vw+\Vwf 
Jn Jn Jn 



\Vv\^ - / iVwP 



2 / IV7 |2 

Jn 
So ||A/'o0||hi/2(s) < 11011^1/2(5) for all /i G if^/^^E) and hence ||iVo|| < 1. D 

It is well known that Landweber iteration, together with a stopping principle for the it- 
erations, is a regularization method for solving equation (3) [EHNOO]. Translating these 
standard results to Kozlov-Maz'ya iteration we obtain the following (which can be de- 
duced, except for the statement concerning optimal convergence rates, from the original 
paper of [KM90]). 



Proposition 2. Let u be a solution of the Cauchy problem (1) with Dirichlet data a and 
Neumann data r. Suppose (r^) are approximations ofr such that \\t — t^\\h-t-/2(s) < ^- ^^^ 
ipo E H^^^'^{B) and let ^^ be the first Kozlov-Maz'ya iterate for the data (a, r^) starting 
from ipQ such that 

where X > 1 is a fixed constant. Then 

Moreover, the rate of convergence is order optimal. That is, if dnu\B = {NqNqY{^) for 
some n > and some ^ E H^^^'^{B), then 

\\dnu\B - ^iWn^^n^B) < c5^'^/(^^+^)eV(2.+i) 
where E = \ |^| |h-i/2(b) (^>^d where c is constant independent of the sequence (r^). 

The previous result assumes that the Dirichlet data a is known exactly. If a is only ap- 
proximately known, this corresponds to error in the operator A^. It is straightforward to 
transform this error into increased uncertainty in the right-hand side of equation (4). To do 
this, we define Fdn ■ H^'^{S) -^ H-^'^{S) by ^^^(7) = dnz\s where 

-Az = \nVL 
z = 7 on 5" 
dnZ = onB. 

A simple computation (left to the reader) shows the following. 



Lemma 5. Suppose that a^ E H^^'^{S) satisfies ||o" — ae||j|^i/2(5) < e. Let N*" be the 
corresponding operator in equation (3), and let r" = FoN^a'' — a). Then for all ip G 

As a consequence, the corresponding termination condition for Kozlov-Maz 'ya iteration 
when there is error in both r and a should be adjusted to 

Ik*^ - V'n||ii--i/2(s) < A(5+ 1 1 r' 1 1^1^-1/2(5)) < K^+ ||-^DAf||e)- 

It is worth remarking that this is an inconvenient criterion to work with in practice: it 
requires both an estimate for operator norm of the Dirichlet-to-Neumann map Fdn as well 
as the size of the error in r in the space H~^/'^{S), which has a rather abstract norm. In 
fact, in many applications (including the work in [MTAS08] that motivates this paper) the 
Neumann data r is known exactly (e.g. r = 0) but there is error in the Dirichlet data a. 
Hence we now consider the Dirichlet version of operator equation (3). 



4 Landweber iteration (Dirichlet version) 

In the previous section we recast Kozlov-Maz'ya iteration with operator KM as a form of 
Landweber iteration by considering a map from Neumann data on B to Neumann data on 
S with fixed Dirichlet data on 5*. The dual formulation obtained by swapping the roles of 
Dirichlet and Neumann data corresponds to Kozlov-Maz'ya iteration with operator MK, 
but posing it requires a little care. The natural operator equation to consider is 

D{ct>) = a (10) 

where D[(j)) ='D{(j))\s. Defining Dq using Vq in place of V, we rewrite this equation as 

Do(0) = a-D(O). 

We would like to consider Dq : H^^'^{B) — )■ H^/'^{S), but the challenge is to find inner 
products on these spaces such that the resulting adjoint D^ leads to a lemma analogous to 
Lemma 2. Unfortunately, this is not true for the norm defined in Section 1.1.3, and it is not 
clear how to adjust it to remedy this situation. 

To circumvent these difficulties, pick any fixed a E H^^'^{dVL) such that (t\s = a and 
lei (f) = a\B. For example, one can obtain such a by applying MK to any element of 
H^^'^(B). Writing = + r^ for some i] E H^^'^{B), equation (2) implies equation (10) 
can be rewritten 

Doir]) = ^-D{^). (11) 

The gain here, as proved in the following lemma, is that the right-hand side of equation 
(11) belongs to H^q^S), not just H^^'^^S), and that if a solution exists, then t] E hIq{B). 

Lemma 6. 
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. Ifr^ e Hli\B) then D,{r^) E Hli\S). 

• Ifcj) = (p + r] solves equation (10) then r] G //qq (B). 

^ 1/2 

• The distribution a — D{(j)) belongs to H^q (S). 

Proof. The claims all follow from the following observation: if (pi and (j)2 are distributions 
in H^^'^{S) that admit extensions cpi and 02 in H^/'^{dVL) that are equal on B, then 01 — 02 G 
ifi/2(5). Indeed, (0i - 02)|i? = and (0i - 02)|5 = <Pi - 02, so 0i - 02 G i/io'(5). A 
similar result holds interchanging S and B. 

Now suppose Tj G i^Qo (-B). Since and Do{ri) both admit H^^'^ extensions to d^ that are 
equal to r] on B, it follows that /^ol??) = Do{r]) - e Hq^'^^S). 

Suppose (J) = (f) + r] solves equation (10). Since and both admit H^^'^{dQ) extensions 

/v I/O 

that are equal to a on S", it follows that ?7 = — G -ffpo (-^)- 

Finally, recall that is the restriction of an H^^'^^dVL) extension of a to B. So a and -D(0) 
admit H^/'^{dVl) extensions that are equal to on B. Hence a — D{(f)) G Hqq (S). D 

As a consequence of the previous lemma, we can interpret equation (11) as an operator 
equation from hI^^B) to hII\S). For t] G Hq^^B), let 

Lniv)=V + D;{{a-D{4>))-Do{v)). 

1 /*? 1 /'? 1 /"? 

Here we treat Dq as a map from Hqq [B] to Hqq {S) and Z^q as a map from /Jqq (5) to 
i^QQ (-B). Landweber iteration (with relaxation constant a = 1) applied to equation (11) 

1 /9 

corresponds to starting with an initial estimate r]o G H^q (B) and computing subsequent 
iterates 77^+1 = Lj:){rik). We then obtain approximate solutions 0^ = + rfk of equation 
(10). We will show that the iterates 0^ are exactly the iterates produces by Kozlov-Maz'ya 
iteration with then operator MK starting with the initial estimate 0o = + 770- To do this, 
we first compute the adjoint of Dq and then prove analogues of Lemmas 2 and 3. 

Lemma 7. Let 7 G H^q {S), and let q and r be the solutions of the boundary -value prob- 
lems 

-Aq = -Ar = inVt 

q = ■y r = on S (12) 

g = (9„r = —dnq on B. 



Then 



Dl{l)=r\B. 
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Proof. Let (j) G i^oo (-^) ^^'^ 7 ^ -^oo i^) be arbitrary, and let v and w be solutions of the 
following boundary- value problems: 

-At; = -Aw = inf] 
(9„f = ti; = on S* 

V = (j) w = (J) on B. 

The equation for w is well-posed since G i/go ("S") and hence the prescribed boundary 
values lie in H^/'^{dVL)\ a similar remark holds for q in equation (12). 

Notice that v — w and q are harmonic, equal zero on B, and equal Do{(f)) and 7 respectively 
on S*. By the definition of the inner-product on Hqq (S) we conclude that 



-o(^) >n 



{Do{<f)), 7)^1/2(5) = I (V(t; - w), Vg) . 
But 



Vv-Vq= / (dnv) q = 
n Jan 

since g = on S and dnV = on S". On the other hand. 



— / Vw ■ Vg = — / w dnq = w dnT = / Vw ■ Vr, 

Jn JdQ JdQ Jq 

where we have used the fact that w = on S and 9„r = — (9„g on B. Since w and r are 
harmonic, equal zero on S, and equal and r\B respectively on B we have 



\7w ■ Vr = {(j),r\B) „i/2 
n 



K!> (B) 



Combining all of the equalities seen thus far we conclude 

for all G HqI'^{S) . Therefore D* (7) = r | b- □ 

-1 /o 

Lemma 8. For any m ifgo (-B), 

D*Do{(l)) = (P-MKo{(P). 

Proof. Let G i^oo (^) ^'^d let u, v, and w be solutions of the following boundary-value 
problems: 

-Am = -At; = -Aw = infi 

S^M = f = w = onS* 

u = (p V = (j) dnW = dnU on B. 
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Let q = u — V and r = v — w. Then q and r solve 



g = 


-Ar = 




in i7 


q = u s 


r = 




on 5 


q = 


a„r = dnV 


-dnU 


on 5 



Noting that mJ^ = -Do(0) and (9„f — 5„m)|b = — <9nq'|B it follows from Lemma 7 that 

L'*L>o(0) = r|B = (t;-w;)|B. 
But v\b = 4> and w\b = MKo{(f)). Hence 

D*Do(0) = 0-Miro(0). 



D 



Lemma 9. 

Proof. Let m, t>, and w be solutions of the following boundary-value problems: 

~Au = f -Av = f -Aw = f infi 

u = cr dnV = T w = a on S 

U = (j) V = (j) dnW = dnV OU B. 

Notice that v\s = D{(p) and w\b = MK{<j)). Let q = u — v and r = w — u. Then q and r 

solve 

-Ag = -Ar = infi 

g = a - D{4>) r = on S 

q = dnT = —dnq on B. 

Lemma? implies Dl{a — D{(f))) = r\B = {w — u)\b- But w\b = MK{(f)) and u\b = 0. 
Hence 

D;{a-D{4>)) = MK{4>)-4>- 

D 

The proof of the following proposition exactly follows the proof of Proposition 1 using 
Lemmas 8 and 9 in place of Lemmas 2 and 3. We omit the proof. 

Proposition 3. For any r] G //qo i^)' ^^ ^'^^^ 

Consequently, the iterates produced by the (Dirichlet) Landweber method and the (Dirich- 
let) Kozlov-Maz 'ya alternating method are identical. 
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Just as with the Neumann formulation, the operator norm of Dq is bounded above by 1, 
which justifies setting the relaxation constant a = 1 in our definition of Ld- 

Lemma 10. The operatornorm of Dq satisfies \\Do\\ < 1. 

Proof. Let rj G i^oo (-^) ^^'^ ^^^ ^' ^' ^^'^ ^ satisfy the following boundary- value problems: 

-Au = -Af = -Aw = infi 
dnU = V = w = u on S 

u = 7] V = f] w = on B. 



Then DqT] = u\s and 



hU\b) = I IVH' 



ll^°^ll^o^o^V)=i^l^^ 



2 



Notice that m — f is harmonic, equals on B, and equals u on S. Hence u — v = w. 
Moreover, 

/ Vu ■ Vf = / {dnU) V = [dnU) TJ = {dnU) U = \'Vuf 

Jn Jdn Jb Jan Jn 

since dnU = on S* and v = u = i] on B. Hence 

-1 /O 

So ||-Do'7||rTi/2,c.N < |k|Li/2.oN for all r] G i/no (B) and consequently ||-Do|| < 1- □ 

-"00 (.'-'> -"do (-^I 

Standard results for Landweber iteration (interpreted in the language of Kozlov-Maz'ya 
iteration) imply the following analogue of Proposition 2. 

Proposition 4. Let u be a solution of the Cauchy problem (1) with Dirichlet data a and 
Neumann data r. Suppose (a^) are approximations of a such that lla — o""^!! i/2,„-, < 6. 

-"oo \^' 

Let (po E (f) + Hqq (B) (i.e. let 0o admit an extension in H^/'^iVt) that equals a on S) and 
let 0^ be the first Kozlov-Maz'ya iterate for the data {as, r) starting from (po such that 

\\a - 0„|| 1/2, , < X6 
where X > 1 is a fixed constant. Then 

lim||«|5-0^||^V2(^^ = O. 

Moreover, the rate of convergence is order optimal. That is, ifu\B — 4> = {D'I^DqY^S,) for 

1 /9 

some fi > and some ^ G i/go i^)' ^^^" 

where E = \ {(.Wfj-^/^r^) <^nd where c is constant independent of the sequence {cr^). 
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The previous result assumes that r is known exactly. This actually holds in many appli- 
cation of interest where r = represents a stress-free or perfectly insulating boundary 
condition. In particular, it holds in the motivating problem from in [MTAS08]. If r is only 
known approximately, then error in r can be rewritten as expanded error in a in a procedure 
analogous to the one described in Lemma 5. 

A more serious weakness of Proposition 4 is that it assumes I la — cr^l I „i/2,„, — t- which 
morally implies that the values of a at the interface of S and B are known exactly. We 
would prefer to have a theorem treating the case \\o- — as\\H^/2is-j — )■ 0. Nevertheless, 
Proposition 2 has some application in this case as well. 

Suppose cr"^ —7- cr in H^/'^{S), and let u^ and v^ be solutions of the following boundary- 
values problems 

- A M^ = - A t;^ = in n 

u^ = a^ dnV^ = dnu^ on S 

dnU^ = v^ = Q onB 

and let m° and v^ be the solution of these problems with a^ replaced with its true value a. 
Then u = u^ — v^ + w solves the Cauchy problem 



Au = f 


in fi 


dnU = T 


onS 


u = a^ 


onS. 


roblem 




Aw = f 


in f] 


dnW = T 


on 5 


W = V^ 


on^. 



if and only if w solves the Cauchy problem 

/ inf] 

(13) 



Noting that v^ G HqI^^{S), and that v^\s -^ v'^\s in H^''^{S), Proposition 2 can be applied 
to the Cauchy problems (13). The stopping criterion then involves the operator norm of the 
map taking a^ G H^''^{S) iov^\se H^i^S). 



5 Conjugate gradient alternative 

Since the Kozlov-Maz'ya alternating method is simply a form of the Landweber method, it 
becomes clear how it might be effectively accelerated. One standard, attractive choice is to 
use the conjugate gradient method. This strategy, together with the Morozov discrepancy 
principle, provides a fast, order-optimal regularization scheme (see, eg., [Ha95]). 

For definiteness we treat the Dirichlet case and consider the normal equation 
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The conjugate gradient algorithm for this problem then reads 

1 ro = a - D{(l) + r]o); 

2 go = D^ro; 

3 d = qo; 

4 k = 0; 

5 while true do 

6 s = Dod; 

8 Vk+1 = Vk + ad; 

9 rfc+i = Tfc - as; 
10 Qk+i = D*rk+i; 

12 d = Qk+i + Pd; 

13 k = k + 1; 

14 end 

Algorithm 1: Conjugate gradient version of Dirichlet Landweber approach 

Using the discrepancy principle, the main loop is terminated when Ikfcll r7l/2.c^ is suffi 



00 



ciently small, and the regularized solution of the Cauchy problem is then V{(j) + rjk). 

The computation of a and /3 requires computation of several norms in i^QQ (B) and H^q (S) , 
each of which would appear to require the solution of a boundary-value problem. We cir- 
cumvent this difficulty by representing each element 7 of Hqq {B) (i.e. qt and d in the 
algorithm) by a harmonic function that is equal to 7 on i? and equal to on 5; the norm 
is then easy to compute according to the definition in Section 1.1.3. A similar principle 

1 /2 

applies to the variables r^ and s in Hqq [S) which are represented by harmonic functions 
that equal zero on B. For this convention to be effective, we need to be able to compute the 
action of Dq and Dq, which the following lemmas show is remarkably easy. 

Lemma 11. Suppose w is harmonic and equals zero on S. Then 

z = Vo{w\b) - w 
is harmonic, equals Dq{w\b) on S and equals zero on B. 

Proof. Since 2; is a difference of harmonic functions it is harmonic. Since w\s = we have 

z\s = ^0(1013)18 = Do{w\b) 

by definition of Dq in terms of "Dq- On the other hand, 150(^15)15 = w^ls, by the definition 
ofpQ. So^Ib = 0. D 

Lemma 12. Suppose w is harmonic and equals zero on B. Then 

Z = -^fo{^nW\s) (14) 

is harmonic, equals Dq{w\s) on B and equals zero on S. 
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Proof. That z is harmonic and equals zero on S follows immediately from the definition of 
A/q. On the other hand, inspecting Lemma 7 with w playing the role of q and z playing the 
role of r in equations (12) we see that z\b = Dq{w\s)- □ 

It is worth remarking that z in equation (14) satisfies the weak formulation that z\s = Q 
and 



V^ ■ Vx = / Vw ■ Vx 

for all test functions x that equal zero on S. Hence the exterior derivative dnW need not be 
explicitly found when solving for z. 

^ ^ 1/2 

Starting with an initial value + r^o G + Hqq (B), let uq be the solution of 

— Auq = f in i7 

Uq = (7 on S 

Uq = (f) + r]o on B. 

Then Algorithm 1, can be rewritten as follows. 

1 ro = uo -V{uo\b); 

2 go = -J^oidnTolB) ', 

3 d = qo; 

4 k = 0- 

5 while true do 

6 s = Vo{d\B) - d; 

8 Uk+i =Uk + ad; 

9 Tfc+i = Tfc - as; 

10 g^+i = -Wdnrk+i\s); 

11 /3 = LlVgfc+i|YL|Vgfc|'; 

12 d = qk+i + Pd; 

13 k = k + 1; 

14 end 

Algorithm 2: Simplified Dirichlet conjugate gradient approach 

When the main loop is terminated (e.g. using the discrepancy principle), the regularized 
solution of the Cauchy problem is then V(uk\B)- Each iteration of the loop requires solving 
exactly two boundary-value problems (one for Vq and one for Mq), just as for Kozlov- 
Maz'ya iteration. In Section 7 we demonstrate how the number of iterations needed for 
the conjugate gradient algorithm can be substantially less than standard Kozlov-Maz'ya 
iteration. 
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6 Variations of the conjugate gradient approach 

We have worked with solving the equation 

D(0) = a 

where D : H^/'^{B) — )• H^/'^{S). By changing the source or target spaces to be L^ spaces, 
one obtains three alternative possibilities for the conjugate gradient algorithm. 

. [L\B) ^ L\S)] 

This variation was treated in [HLOO] . Because the choice of domain has lower reg- 
ularity than H^/"^, one expects the reconstructed solutions to exhibit lower regularity 
than Kozlov-Maz'ya iteration. Indeed, one step of the algorithm involves a boundary 
condition of the form 

u\b = dnV\B, (15) 

where f is a previously computed harmonic function. This step has the effect of 
lowering the regularity of u and is perhaps responsible for oscillations observed in 
[HLOO] Figures 2, 4, and 6. In particular, we performed a reconstruction of [HLOO] 
Figure 2 using the H^^'^ — ;■ H^^'^ method (as well as the H^^'^ — ^ L^ method de- 
scribed below) and these oscillations are absent. 

• [L'^{B) —^ //^/^(S)] This approach appears in [Kn04], although it is not presented 
as such. That paper considers the functional 



where v is the solution of 



G{(P) = / \V{v-w) 
Jn 

-Av = f in fi 

V = (7 on S 

V = 6 on B. 



|2 



and w is the solution of 



— Aw = f in i7 
dnW = T on S 
w = (p on B. 

Noting that (v—w) is harmonic and equal to zero on B, we see that G can be rewritten 
as 

G(0) = \\iv-w)\s\\l,^,^^^ = \W-w\s\\l.f.^,y 

which is the functional being minimized by the Dirichlet Landweber procedure pre- 
sented in Section 4. However, [Kn04] initially uses the L'^{B) gradient of G. Since 
G is not defined on all of L^ (we do not expect solutions with L^ boundary data to 
lie in H^{n)), the L^ gradient leads to a loss of regularity, and the algorithm has 
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a step similar to equation (15). This trouble is ameliorated in [Kn04] by introduc- 
ing a smoothing step, effectively recasting the domain as a subspace of H^{B) and 
factoring the map through L'^{B). Since the norm used for H^{B) involves a PDE 
defined only on the domain boundary, the additional step adds some complication to 
the algorithm when compared to Kozlov-Maz'ya iteration. 

{H^'^{B) -^ L2(5)] 

This combination does not appear to have been previously addressed in the literature, 
and has some potential interest. Since the domain is H^^'^{B), we expect a recon- 
struction with higher regularity than the method of [HLOO]. And since the range is 
L^(5'), the associated stopping principle will involve LF' error estimates, which are 
much easier to obtain than the rather abstract H^/"^ error estimates. 

As in Section 4, we consider the operator equation 

{ioD){r,)=a (16) 

where D : H^/'^{B) -^ H^/\S) and where i : H^/\S) -^ L\S) is the natural 
embedding. Note that we work with H^^^(B) rather than the more awkward space 
i^oo {B). Equation (16) can be rewritten 

and to apply the Landweber or conjugate gradient methods we need to be able to 
compute (l o Dq)*. 

Proposition 5. Let 7 G L^{S), and let w be the solution of 



- 


-Ait; = 


in VL 




dnW = 7 


on S 




dnW= ^1 


/ 7 on B 


1 

B\ 


j W= J. 

J B Js 





(17) 



Then {ioDo)*{-f) = w\b. 

Proof. We first note that there is a solution w E H^{ri) of system (17). Indeed, the 
Neumann data (which belongs to L'^{dVL) C H^^^'^{dVL)) satisfies the compatibil- 
ity condition j^^ dnW = 0. Hence the PDE admits a solution in H^{n) determined 
uniquely up to a constant. The final equation then determines the value of the con- 
stant. 

Now let G H^^'^{B) and 7 G L'^{S) be arbitrary. Let w be the solution of system 
(17) and let v and q solve 



At; = 


- A g = in fi 


dnV = 


dnq = on S 


V = (f) 


q = w\b on B 
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Then -Do(0) ="^15 ^^'^ 



((^°^o)(0),7)l2(5) = / ^7 



/ V dnW + I V 

'dn Jb . 



\B\ 



1 

Is J 



Vv-Vw + 



Now 



Moreover, 



7 



1 

LBI 



w 



B 



1 

LBI 



1 

\b\ 

B 



IB J 



f 7 
Us J 



Vm ■ Vw = / dnU w = dnU Q = Vm ■ Vw 
/n J an JdQ Jn 

since 9„m = on S and q = w on B. Combining all these equations we conclude 



((ioDo)(0),7)L2(5) 



Hence (io Do)* (7) = qIb- 



Vv-Vq + 

{v\B,q\B)m/2(^B) 



1 

LBI 



B 



1 

LBI 



q 



B 



D 



The conjugate gradient algorithm for this problem starts with an initial value 0o £ 
H^^'^(B) and a corresponding uq = V^cpo). Given 7 G L'^{S), we define Wo(7) = 
w|b where w is the solution of system (17). We then obtain an analogue of Algorithm 
2 by tracking elements of H^^'^(B) as harmonic functions with zero Neumann data 
on^. 

1 ro = a-Uo\s; 

2 go = Vo{Wo{ro)) ; 

3 d = qo; 

4 k = 0; 

5 while true do 

a = /f^|Vgfc|Y /sMIsI^; 
Uk+i = Uk + ad; 
rk+i = rfc - ad\s; 
Qk+i = VoiWoin)); 

d = qk+i + f3d; 
k = k + 1: 



11 
12 
13 end 



Algorithm 3: H^/'^{B) — )■ L^(S') conjugate gradient approach 
Upon exit, the regularized solution of the Cauchy problem is simply Uk- 
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Figure 1: Left: computational domain. Right: True values of a = u\s and (p = u\b. 

Kozlov-Maz'ya 

if 1/2 ^ H'/2 




iteration 



Figure 2: History of the L^ discrepancies for the p = 0.1 reconstruction. 

7 Numerical Results 

Let n be the domain in the plane bounded above by the x-axis for — 1 < a; < 1 and bounded 
below by a parabola passing through the points (—1, 0), (1, 0), and (0, —d) (Figure 1, left). 
We take S to be the region on the x-axis with — 1 < x < 1 and B to be the portion of the 
boundary with y < 0. The Cauchy problem to solve is 

— A M = /o in i7 

u = a on S 

dnU = on S* 

where /o is a constant. This is a model for a glaciological inverse problem where Vt is the 
cross-section of a glacier and u represents the component of ice velocity orthogonal to the 
cross-section. The homogeneous Neumann condition arises as a consequence of a zero- 
stress hypothesis at the ice surface S, and surface velocity measurements are represented 
by a. 
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Figure 3: Reconstruction of u on B. Left: p = I. Right: p = 0.1. Kozlov-Mazya (black, 
dotted), H^^^ — )■ if ^/^ (red, dashed), H^^^ — > L^ (green, solid), true value (blue, solid). 
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Figure 4: Detail of Figure 3, p = 0.1, near the interface of S and B. 



We consider synthetic data a obtained by numerically solving the problem 



Au = fo in r2 
onS 
u = 6 on B 



dnU 



(18) 



where is a prescribed function; then a = u\s. We used a bump-function 

(f){x, y) = Mo exp ( -- x^/s^ 

where ?io and s are constants. Figure 1 (right) illustrates values of ?i on B and S. 

All computations were done using the finite element method (and in particular using the 
FEniCS [LWIO] framework). In the following we used depth d = 1/2, forcing term /o = 8, 
peak basal speed mq = 1/2^ and standard deviation s = 1/3. The forcing term was selected 
so that the peak value Wmax of the solution u of system (18) was approximately 1. 

The surface measurements a were perturbed by spatially uncorrected gaussian noise with 
standard deviation Wmaxioo^ where p is a constant describing 'percent-noise'. We then 
applied standard Kozlov-Maz'ya iteration, the H^^"^ — )■ H^^"^ conjugate gradient method 
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(Algorithm 2), and the iJ^/^ — )• L^ conjugate gradient method (Algorithm 3) to solve the 
inverse problem. The starting estimate was 0o = for all reconstructions. 

Estimating error in the H^^'^{S) norm for use in the discrepancy principle is impractical, 
whereas we have good estimates for the error in the L'^{S) norm. For all three algorithms, 
we therefore terminated early based on the L^ discrepancy. This is formally justified only 
for the H^^^ — )■ L^ algorithm. Nevertheless, the other two algorithms using this modified 
discrepancy principle appeared to remain stable. It would be interesting to have a formal 
proof of this observation. 

Figure 2 shows a graph of the L^ discrepancies 

\\u-Uk\\LHS) 

for the p = 0.1 reconstruction for each of the three algorithms. As would be expected, 
the conjugate gradient based algorithms required substantially fewer iterations to reach 
their target discrepancies. This improvement was less substantial for the larger (and more 
physically relevant) value of j9 = 1 where the conjugate gradient algorithms each used six 
iterations and Kozlov-Maz'ya iteration used nine iterations. 

All three algorithms gave qualitatively similar reconstructions for p = 1 and p = O.l 
(Figure 3). The algorithms using H^^'^^S) for the surface norm tended to have stronger 
oscillations near the interface of S and B. Figure 4 shows a detail of these oscillations near 
the boundary point (1, 0) in the p = 0.1 case. Although present in all three reconstructions, 
the oscillations were more damped by the H^^^ —^ L^ algorithm. 

Of the three algorithms used, the H^^^ — )■ L^ algorithm had a number of slight advantages. 
Jn addition to having the speed of the conjugate gradient algorithm, it has a provably rigor- 
ous stopping principle for an easily obtained error estimate, and it generated subtly better 
reconstructions near the interface of S and B. 



Appendix 

We recall here standard facts about Sobolov spaces. Unless otherwise noted, we use the 
definitions and notation of the careful exposition in [McOO]. 

Let i7 C M" be a connected, relatively compact open set with a Lipschitz boundary dD.. 
Suppose S and B are open subsets of dQ sharing a common boundary 11 and that dil = 
5 U 5 U n is a Lipschitz dissection as defined in [McOO]. 

For s E [—1, 1], we define H'^(di^) as in [McOO]. In particular, the trace map from if ^(^2) 
to H^/'^{dVL) is a continuous surjection that admits a continuous right inverse. 

The space H^{S) is then defined as restrictions to S of distributions in H^{dVL), and is 
given the quotient norm, so 

\\u\\h-{S) = mi{\\U\\Hs(dn) ■ U\s = u}. 
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We define Hqq(S) to be the closure in H^(dQ,) of the Lipschitz functions with compact 
support in S* . In particular, Hqq(S) is a closed subspace of H'^{di^), though we will 
commonly identify such functions as elements of H'^{S). Because of the regularity of the 
common boundary 11, 

H^oiS) = {ue H%S) : u\b = 0}. 

Spaces of distributions on B are defined similarly. 
The bilinear form on Lip(9f2) x Lip(9i7) given by 



{u, v) = I uv 

extends continuously to a bilinear form on H^'^(dri) x H'^(dfi), and the map u i-> {u, ■) 
gives an isomorphism between H'^^^dVl) and H^{dri)*. Again, because of the regularity of 
the common boundary 11, there is an isomorphism l between Hqq^{S) and {H'^{S))* defined 
as follows: for u G Hq^^S) and v G H%S), 

iLu)iv) = {U,V) 

where V E H^{dri is any function with V\s = v. By reflexivity we also have an iso- 
morphsim between H^S) and Hoq{S)*. Note that [McOO] denotes ifg'o by H'. 



Mixed Boundary Value Problems 

We wish to solve 

— Au = f in n 

u = (f) onS (19) 

dnU = T on B. 

Let H^{n, S) = {u E H'^i^l) : u\s = 0}. This is a closed subspace of H^{n) (being the 
kernel of restriction to the boundary followed by restriction to S). For u,v E H^{ri, S) let 



A{u,v) = / Vn- Vv. 
Jn 



This is evidently a continuous bilinear form. Moreover, it is strongly coercive so long 
as S* is nonempty (since we have assumed additionally that S is open). The argument is 
completely analogous to the corresponding one for the well-known case where S = dVl. 

Now suppose / E {H^i^))*, (j) E H^^^S), and r E R-^/^B). A weak solution of 
equation (19) is a function u E H^{i^) such that u = (p in S and such that 



A{u,v) = fv+ TV (20) 

Jn Jb 

for all V E H^{ri,Ti). The integrals on the right-hand side of equation (20) are to be 
interpreted as shorthand for the application of linear functional. In particular, the integral 
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over B makes sense for if v G H^{Q, S) then t>|5 = and hence v\gQ G H^^q (B) = 
H-^I\BY. 

Pick w E H^iTl) such that w\s = 4>', indeed since the trace and restriction operators have 
continuous right inverses, we can pick a w such that its norm in H^iVt) is controlled by the 
norm of in H^/'^{S). Then u = w + h solves equation (20) if and only if h E H^(Vl, S) 
and 



\/h-Vv = / {fv - Vw ■ Vv) + TV (21) 

n Jn Jb 

for all V E H^{fi, Fi). The right-hand side of equation (21) is a continuous linear functional 
on V G H^{ri, Fi). So the the Lax-Milgram theorem implies that there is a unique solution 
h of equation (21) and hence a unique solution u of equation (20). Moreover, the Lax- 
Milgram theorem (along with the aforementioned control on the size of | |w| |j:/i(q)) implies 
there are constants ci and C2 such that 

\\u\\m{n) < ci{\\w\\H-i{Q) + \\h\\m{n,s)) 

< C2i\\f\\{HHn,S)r + |l0|l/fi/2(S) + lklU-V2(B))- 

Boundary Neumann Data 

Suppose u E H^{Q,) and — Au = f where / G (iJ^(f2))*. Then we define dnu\dn ^ 
H-^/\dn) by 

/ (9„m)0 = [vu-Vv~ [ fv (23) 

JdQ Jn Jn 

where v is any element of H^(n) that equals (p on dVl. See, e.g., [McOO] for a proof that 
this is well defined. We repeatedly use the following lemma in various guises (and with 
little comment) throughout this paper. 

Lemma 13. Suppose v and w belong to H^{Q), — Aw = f E {H^{Q))*, v\s = 0, and 
dnwls = 0. Then 

/ V dnW = 0. 

JdU 

-I /o 

Proof. Since v\s = we have v\dn ^ Hqq {B). So there is a sequence (0^) of Lipschitz 
functions on dVt with supports contained in B converging to v in H^/'^{dVL). By definition 

Js 'Pk dnW = fg^ 4>k dnW. Siucc dnw\s = 0, fg 0fc d„w = 0. So Jq^ 0fc dnW = for each 
k. Since dnW E H-^/^{dQ) = (H^^^idQ))* and since cf>k ^ v in H^l^{dn) we conclude 
that Jg^ V dnW = 0. n 



Equivalence of Norms 

We sketch the proofs here that the norms described in Section 1 are equivalent to the stan- 
dard norms on those spaces. In the following we use three-barred norms ||| ■ ||| to denote 
those from Section 1 and reserve two-barred norms for their standard definitions. 

25 



Let ip e H ^/^(S*) and let u be the solution of 

- A M = in fi 

dnU = ijj on S 

u = on B. 

So lll'^lll_ffi/2(s) = ||V'u||i2(Q). We see from equation (23) that 

\\drM\H--^/'\dn) < l|Vn||i2. 
By definition, 

\\dnU\\H-i/2^s) < II'9„w||h-i/2(0) 

and hence ||'?/'||j|^-i/2(5) < ||Vn||L2(Q) < ||u||i:/i(Q). On the other hand, equation (22) 
implies 

\\u\\m{n) < C2||V'||h-i/2(5) 
for some constant C2 independent of V^. Hence 



which establishes the desired equivalence. 

-1 /o 

Let G Hqq (S) and let u be the solution of 

- A M = in fi 
u = (j) on (9^2 

so lll0lll//l/2.5,^ = ||Vm||^2- From inequality 22, the continuity of the trace map if ^(r2) — )■ 

H^/'^{dVl), and the definition of the norm on i^QQ {S) we see that there are constants Cq and 
C2 such that 

Co||0||hi/2(5) < ll-^ll/ZHf^) ^ C2||0||^i/2(5) 

But for functions in H^ vanishing on B it is well known that ||u||Hi(n) is equivalent to 



Similar arguments work for the equivalence of norms in H^^'^^S). The only new ingredient 
is the fact that the norm | |u| |/^i(q) is equivalent to 



|Vm|^ + 



1 f ^^'^'^ 
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